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ABSTRACT 


Wave blocking in river inlets is examined using the NHWAVE (Non-Hydrostatic Wave) 
model under development. Blocking flows at river inlets are a significant hazard to 
navigation. Refractive and shoaling effects contribute to the enhancement of wave field 
energy, causing instabilities and breaking, resulting in energy dissipation and transfer at 
the blocking point. The non-linearity of wave-current interactions and wave breaking 
makes the dynamics of blocking flows difficult to model. Current efforts to use wave- 
averaged models are insufficient to describe the complex dynamics that occur within one 
or two wavelengths of a blocking point. NHWAVE uses the non-hydrostatic, 
incompressible Navier-Stokes equations to model fully dispersive wave processes in the 
time domain. Monochromatic wave cases are explored and compared with lab 
experiments of energy dissipation due to wave breaking under conditions of strong 
opposing current, conducted in 2002 at the University of Delaware by A. Chawla and 
T. J. Kirby. The model was initially unable to resolve the boundary conditions necessary 
to model wave blocking in a tank domain. However, developments to the numerical 
scheme in NHWAVE have advanced its capability in this regard. Due to the difficulties 
of modeling the dynamics of wave blocking and the boundary conditions in a wave tank, 
only preliminary results were obtained. NHWAVE needs further development; it shows 
promise, however, to be able to predict wave reflection, blocking, and dissipation on a 


strong opposing current. 
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I. INTRODUCTION 


Wave fields propagating into regions of significant opposing flow, such as at the 
mouth of a tidal inlet, can be a significant hazard to navigation for mariners, affect the 
design of coastal structures, and alter the transport of sediment in the near shore region. 
As a wave field propagates into an opposing current, the group velocity reduces leading 
to an increase in wave height and decrease in wavelength. These waves can break causing 
energy dissipation, or the group velocity can go to zero in which case the wave field is 
blocked. If the current is not strong enough and the wave amplitude is relatively small, 
the waves can also be reflected from the blocking point. Energy transfers in blocking 
flow cause the wave environment to become very rough and difficulties in navigating and 


transiting these inlet areas are common. 


A. WAVE BLOCKING THEORY 


The basis of previous research on small-scale wave-current interactions is the 
Doppler-shifted dispersion relation for waves moving on uniform currents (Perigrine 
1976). In the case of a wave field opposing a tidal current, the time and length scales for 
the tidal current flow are much larger than the wave period and wavelength. This allows 
kinematic properties of plane waves to be applied on uniform currents. By applying a 
depth uniform current U and using a reference frame that is moving with the current, the 
equations and solutions for water-wave movement are the same as for no current 
(Jonsson 1990). Previous study of waves moving on a vertically sheared but weak current 
showed that by using the small parameter e=U/c, where c is the wave phase speed, 


through O(e) 
(o@—kU cosa) = gk tanh(kh) = 0° (1) 


where w and k are the fixed reference frame radian frequency and wave number, and o is 
the intrinsic frequency defined by linear theory (Kirby and Chen 1989). The still water 
depth is 4 and a is the angle between wave propagation and the underlying current. This 


is identical to the dispersion relation for arbitrary water depth using the intrinsic 


frequency @-keU. Defining as wave frequency in a reference frame moving with the 


current, the wave phase speed in both frames of reference is related by 


O_o 
—=— +|U|cosa (2) 
[AL Ay 
or 
a-k-U=o (3) 


where wave number, k, remains unchanged in the two reference frames. Application of o 


is dependent on the wave theory used. 


The kinematic wave blocking condition is described using the Doppler-shifted 
dispersion relation for waves and is obtained by assuming depth uniform currents, 
Equation (3), along with an expression for . Differentiating with respect to k and setting 
the group velocity to zero gives 


Oo 
Ty 4 
Ak (4) 


where 


U=|U|jcosa (5) 


Due to the cosa term, blocking only occurs for a between 90° and 270° (Chawla and 


Kirby 2002). 


The graphical solution of the Doppler-shifted dispersion relation, Equation (3), 


shows the limits of different wave-blocking scenarios to be addressed in this paper, 
(Figure 1). Using linear wave theory, Equation (1) becomes o =,/gktanhkh , and is 


plotted as a blue curve for the case of currents directly opposing the wave field, a=180°. 
In Figure (1a), Equation (3) is computed with a constant current speed, U, for different o 
values and plotted as a green, red, and cyan line (@;, @2, and 3). For increasing 
frequency and a constant current, U, the solution evolves from wave-current interactions, 
to a blocked wave field. The @, line is a solution for waves below the blocking frequency 
for the given current. These waves will have interaction with the current, but will not be 


blocked. Where a—k-U=o becomes tangent to Equation (1) at point A, defines the 


blocking solution. This Doppler-shifted solution represents the frequency at which the 
given current will block the wave field. Waves above this curve have no solution; this is 


indicated by 3 in which all waves will be blocked. 


The blocking scenario is categorized in Figure (1b) where the dispersion relation 
is solved at a particular frequency, @2. The horizontal line represents a wave only 
scenario with U=0. Solutions below this line, in the dotted area, will have no current 
interactions. This is the case for solutions from point A to point B. From point B to C, in 
the presence of an opposing current, the wave number will increase and the wavelength 
will become shorter. Here, the solutions are shown in the vertically lined section. As k 
increases further and the group velocity decreases, the complete blocking scenario is 
reached and is represented by point C. There is no solution for waves interacting with 
stronger currents above the @2 line and therefore the solution never increases over this 
line for a blocking current. For current smaller than the wave blocking current but with 
higher k, there is a solution in the wiggly line section. In this case, while the waves are 
moving opposite the flow, the wave energy travels with the current, opposite the wave 


motion. The solution represents reflection and partial wave blocking (Chawla, 1999). 





Figure 1. | Graphical solution to the dispersion relation. a) Solution for constant 
current and different frequencies @ 1, ®2, and 3. b) Solution for increasing 
current, dots=wave only, lines=wave-current interaction, 
wiggly lines=wave-current interactions with some reflection 


B. BASIS FOR STUDY 


The Office of Naval Research has directed research initiatives to study the 
complex wave-current dynamics in the nearshore region. As the Navy moves further 
inshore and upstream, these conditions at river inlets become more important to 


understand. This effort is focused on developing a model that can handle the complex 


dynamics and boundary conditions involved in modeling a tank domain with both waves 


and currents so that future application to the field is possible. 


Most previous studies and theoretical analyses on wave blocking have been 
limited to linear wave theory. Due to the sharp steepening of waves just prior to the 
blocking point, this approach is only valid for small amplitude waves. In locations such 
as a river inlet, where wave heights can be great, an approach that can model the non- 


linear effects is needed. 


1. Motivation 


Wave-current interactions are an important consideration in coastal 
oceanography. Non-linear interactions can cause significant effects to mass transport, 
large amplitude wave crests can cause damage to coastal structures, and the turbulence 
induced by the kinematics can be dangerous for maritime vessels. Modeling these effects 
to predict the dynamic conditions in areas of opposing flow and wave-blocking is of great 
importance and is difficult. Despite recent studies in non-linear wave-current interactions, 
there is still a limited understanding of how energy is transferred or dissipated at wave 


blocking locations. 


Field measurements at river inlets, where a strong ebb tidal jet meets an opposing 
wave field inducing significant wave-current interaction and wave-blocking, are 
becoming more common with improvements in field observation techniques and 
instruments. Simulating these conditions in a model poses issues of scaling and boundary 
conditions that are not present in field data. Simulating a tank condition in a model 
domain further removes the dynamics from the environmental truth but is necessary to 
study individual wave components and improve understanding of energy transfers 


through blocking regions. 


Understanding wave blocking with NHWAVE allows for potentially better 
parameterization with more robust wave-averaged models that run faster. Currently those 
models under predict wave heights and breaking by about one half. This is due to the 
non-linear effects occurring at distances on the order of a wavelength, this causes wave 


solutions in wave-averaged models to fail. 


A model of the wave-current interactions at an inlet can potentially describe the 
fully non-linear processes in wave blocking. A number of difficulties are associated with 
numerical implementations in modeling, which will be described in this research. This 
study will review some of the contemporary attempts to model and predict the physics of 


wave blocking and discuss the development of one model in particular. 


2; Objective 


In the Chawla and Kirby (2002) experimental study, hereafter referred to as 
CK02, monochromatic and random wave fields were created with an opposing current in 
a wave flume. A conservation of wave action numerical model was used to describe 
waves on an opposing current. A simple probabilistic bore model was used to describe 
dissipation by wave breaking to replicate the lab data. While good correlations with 
laboratory data were obtained in this simplified model, extension to a more complex 
model that can handle dynamics in the wave-blocking region in the time domain is 


attempted here. 


The objective of this research is to develop the NHWAVE model (described in 
future sections) to predict wave-current interactions and blocking compared with 
extensive laboratory data described in CKO02. Of particular interest is the amplitude 
envelope of monochromatic waves in cases varying from complete blocking and 
reflection, to total wave transmission through the blocking point. Energy transfers and 
dissipation in larger amplitude breaking wave cases will also be explored in the model. In 
the following, data from the tank experiment is used to develop a model in a refractive, 
monochromatic wave setting. NHWAVE has been modified to characterize wave-current 
interactions in a new and innovative way, showing its potential in this modeling 
environment, and adding relevance to more general modeling of wave blocking at river 


inlets. 


HW. THEORETICAL BACKGROUND 


A. INTRODUCTION 


Previous studies of wave-current interactions have been limited to linear wave 
theory. However, when the opposing current is strong enough, waves become blocked, 
and linear theory breaks down at the blocking point. Both refractive and shoaling effects 
contribute to the enhancement of energy, causing instabilities and breaking that result in 
energy dissipation and transfer at the blocking point. The non-linearity of wave-current 


interactions makes the dynamics of blocking flows difficult to model. 


B. PREVIOUS MODEL APPROACHES 
f. Boussinesq Models 


Non-linear Boussinesq-type wave models have been used to describe the 
refraction, diffraction, shoaling, and breaking of wave fields in the coastal region. The 
original equations by Boussinesq (1872), describes water waves in incompressible, 


inviscid fluid and non-rotating flow and are based on combining the horizontal 
Say ; sees a 
momentum flux and the continuity equations and retaining terms Lacing) =O(1), the 


Ursell Number, where a is amplitude, 


On On (3 1.,8n) 

a a aela an 3" ae)” 
1 is the free surface elevation, g is the gravitational acceleration, and ¢ and x are time and 
horizontal location. Models based on this theory with improved non-linearity and 
dispersion characteristics by including higher order terms have been well tested and are 
efficient models of wave propagation. The non-linear aspects of the improved models 
allow the amplitude dispersion and energy transfer effects in wave blocking scenarios to 
be well represented. Conversely, for short waves, wave dispersion properties are not well 


predicted and this deviation only increases with increasing wave number. The governing 


7 


equations do not include dissipation due to wave breaking, and thus, become invalid 
under any breaking conditions. Attempts to include turbulent mixing and dissipation 
caused by breaking in these models have used artificial eddy viscosity terms (Heitner and 
Housner 1970, Tau 1983, Zelt 1991, and Kennedy et al. 2000 amongst others). This 
approach is useful for weakly non-linear waves, with longer wavelengths. However, this 
is not the case in wave-blocking scenarios, and therefore these methods are not sufficient 
to model wave-breaking on an opposing current and the turbulence that occurs at wave- 


blocking points. 


2. Non-linear Shallow Water Equations (NLSW) 


The non-linear shallow water equations assume a hydrostatic pressure distribution 


that is violated at the point of wave-breaking. In basic form 





di de de 
t IX IX (7) 
du dx dD 
+u—+g—=0 
dt du dx 


where uw is the horizontal current velocity and D=7+/h is the water depth. 


Traditionally, the effect of non-hydrostatic pressure can be modeled using the 
Bousinesque-type approximation and adding higher order terms to the NLSW equations. 
The NLSW equations are then extended with the addition of the vertical motion equation, 
a deviation from hydrostatic pressure, making them equivalent to the incompressible 
Navier-Stokes equation. In this approach, as with previous modeling techniques, there is 
difficulty in simulating the free surface that inhibits accurate prediction of wave-current 


interactions at blocking points. 


3. Renolds Averaged Navier-Stokes Equations (RANS) 


These efforts led to solving the Reynolds averaged Navier-Stokes equations with 
simultaneous computation of the vertical structure and horizontal variations, 
ou, | i=1,2,3 


= rma —po,+ [ou — pu,u (8) 
Be Ox, ge Al Os ax, | ns J=1,2,3 


Ou; 


ae ox, 


where f is representative of external forces and the over-bar indicates time averaging. 
The instantaneous velocity and pressure terms are decomposed into a constant and a 
fluctuating component that are time averaged on the scale of the wave period. In 
both the Boussinesq and NLSW equation methods, uniform pressure and velocity 
distributions are assumed. Also, the RANS method does not employ the simple 


turbulence parameterization as in previous methods. 


4. Free Surface Treatment 


Of importance in all modeling techniques is the treatment of the free surface and 
obeying the balance of pressure forces across the boundary. Methods include the 
Volume-of-Fluid (VOF), Marker and Cell (MAC) and Smoothed Partical Hydrodynamics 
(SPH) methods (Zijlema and Stelling 2011). Drawbacks to these model techniques 
include coarse resolution on the interwave scale, inaccurate velocity computations where 
the pressure boundary condition is not met, and computational inefficiency. Because of 
these inadequacies, it may not be possible to accurately predict the onset of wave 


breaking and energy transfers that occur at wave-blocking river inlets. 


The difficulties in modeling the pressure boundary condition necessitate non- 
hydrostatic approximations at the free surface to achieve the pressure balance at the 
boundary. Stelling and Zijlema (2003) used a Keller-box staggered grid, allowing 
pressures to be defined at the vertical cell faces rather than cell centers. This pressure 
boundary simplification and staggered grid in the vertical direction allow for non- 
hydrostatic approximations to be used at the surface, enabling wave propagation, 
turbulence, and solute transport calculations with fewer layers and a better representation 


of linear dispersion effects than previous models. 
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While the above approximations allow for the use of a non-hydrostatic pressure 
surface condition, models of turbulent conditions, such as wave breaking in the surf zone 
and run-up in the swash region, must be adequately modeled for shock propagation 
(Zijlema and Stelling 2008). This is accomplished by applying the Godunov-type 
approach for shock capturing and discontinuous flow. These schemes allow wave 
breaking to be calculated through a numerical approach rather than being prescribed in 


the modeling process (Ma et al. 2012). 
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Ht. MODEL DEVELOPMENT 


A. INTRODUCTION 


Two models were recently developed that use the non-hydrostatic Navier-Stokes 
equations in conservative form to model fully dispersive wave processes, NHWAVE 
(Non-Hydrostatic Wave Model) (Ma et al. 2012) and SWASH (Simulating Waves till 
Shore) (Zijlema et al. 2011). To study the complex wave-current interactions at blocking 
locations, an effort was begun to look at laboratory data and attempt to replicate it in a 


model environment so that future application to field data would be possible. 


Initially, the SWASH model was tested to replicate the tank experiment data. 
SWASH is an open source code model simulating non-hydrostatic, free-surface, 
rotational flows in multiple dimensions. It is useful in predicting transformations of 
surface waves and rapidly varied shallow-water flows in coastal waters. It implements 
nonlinear shallow water equations and generally describes the complex changes due to 
rapidly changing parameters in a shallow water environment (Zijlema et al., 2011). This 
model has been through numerous benchmark tests and handles nearshore dynamics very 
well. The availability of implementation manuals, and vast examples of model use made 
independent work with the model possible. Ideally, this model would be able to handle 
the dynamics of the CK02 experiment. 


Initial test cases with the SWASH model looked promising. The first steps were 
to define a wave tank with a false wall to replicate a river inlet and induce a wave field 
with an opposing current. SWASH has a user-friendly interface and the process was 
relatively simple. Model runs with a monochromatic wave field imposed at the boundary 
were consistent with expected linear theory. Separate tests to determine the current input 
as a mass flux at the boundary were also successful. However, it was not possible to 
include both dynamic inputs in the model. The boundary inputs of the model could 
handle a wave field or a current, but not both imposed over each other at the same 
location. After repeated attempts to correct the problem through the code and different 


inputs, it became apparent that there was no further way forward and help from the 
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developers was needed. Discussion about the boundary condition issue was initiated with 
model developers; however, this was an issue that would not easily be resolved and other 


options needed to be explored. 


Open communication with peers at the University of Delaware in the Coastal 
Engineering Department facilitated discussion of the newly published model NHWAVE 
(Non-Hydrostatic Wave). Being able to consult with the developers of the model on a 
regular basis was of great importance. The model was published in May of 2012 and has 
undergone numerous benchmark tests. Contrary to SWASH, however, its main focus was 
on the prediction of Tsunamis and application to greater ocean depths than SWASH. 
These models differ primarily in their numerical approach to non-linear ocean processes. 
NHWAVE is developed for wave blocking in this research and described in the 


following. 


B. NON-HYDROSTATIC WAVE MODEL (NHWAVE) 


The NHWAVE model uses the o, coordinate system for both surface and terrain 
features as developed by Phillips (1957). The o, coordinate is defined as 





= _(z+h) 


aa; (9) 


where D is the total depth (A+) and z is measured positive upwards from the still water 
level. The coordinate transformation maps the bottom and surface to constant boundaries 
of a strip of unit thickness. 

1. Governing Equations 


Using the conservative Boussinesq equations in general form 


<+V-0(¥)=s (10) 


is the vector for conserved variables and ©(‘P) is the numerical flux vector function. 
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p= (11) 


O= 12 
Dui+ Dvj+ak ie?) 


(Duu + = eb?) + Duvj+uaok 





Duvi+(Dvv+ 5eP)j + vok 





Duwi + Dvwj+ wok 


where (u,v) are velocity components in (x,y) and w is the vertical velocity in the oz 
coordinate system. The source term S, is composed of three components, bottom slope 


S;, pressure gradient S,, and turbulent mixing S,. 

















S,= S,= S.= (13) 
ope (2 4. dp | DS, 
Ox p\dx do dx DS, 
pot ears _ DS 
oy p\dy dody 
? 1 dp 
ae 


here, p is dynamic pressure only. The k-e closure model is used to calculate the turbulent 


diffusion terms S, , S., S,_ . In NHWAVE, the free surface governing equation can be 
written as 


dD a 1 d 1 
+ (PI udo.)+—“(DI vdo.}=0 (14) 
dt dx\ “0 dy\ +0 

The dynamic pressure in the source terms is calculated by solving the Poisson equation in 


(x, v, oz) coordinate system. Further details of the calculation techniques are found in Ma 


et al (2012). 
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2. Numerical Scheme 


Finite-volume and finite difference methods are used for special discretization. It 
is done in two steps, first using a reconstruction technique to compute values at cell 
interfaces, and second using a Reimann solver for numerical fluxes at the interface 
location. A central difference scheme is used for all source terms, S. NHWAVE uses a 
second-order Strong-Stability-Preserving (SSP) Runge-Kutta scheme for non-linear 


discretization 


vw) =" +Ar(-V-0(¥") +5") 


y= sv +[¥ + Ai(-v-0(¥")+5™) | ie) 


Within the above scheme, a two-step projection method splitting the time integration into 
hydrostatic and non-hydrostatic steps is used. An adaptive time step is used following the 
Courant-Friedrichs-Lewy (CFL) criterion. NHWAVE is run parallel using the domain 
decomposition technique. Non-blocking data communication between processors is done 


using the Message-Passing-Interface (MPI) when not in serial mode (Shi et al. 2012). 


Computationally, a staggered grid approach is used defining velocities at cell 
centers and pressure at the vertical facing wall similar to the Keller-Box scheme. This 
allows accurate treatment of the pressure boundary condition at the surface. NHWAVE is 
capable of simulating wave refraction, diffraction, shoaling, breaking, landslide tsunami 
generation, and alongshore current. It can predict surface wave processes using few 
vertical layers with good accuracy as seen in benchmark tests (Tehranirad et al. 2012). 
Wave breaking and associated energy dissipation are also reasonably predicted by the 
model. Additional details of the development and testing of NHWAVE are found in Ma 
et al. (2012). 


C. INITIAL MODEL TESTING 


There are limitations inherent to ocean modeling and numerical code. While 


modeling of the physical processes has become computationally simpler with 
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developments to wave theory, any computer simulation will have difficulty in modeling 


all variables and especially boundary conditions. 


There is a world of difference between the character of the fundamental 
laws, on the one hand, and the nature of the computations required to 
breathe life into them, on the other. Berlinski (1996) 


The current study was conducted through several iterations of the NHWAVE code. Both 
modifications to the source code, compilers, and the executable files were needed to 
model the wave-current interactions in a tank. NHWAVE continues to be under 


development and shows promise for future modeling efforts. 


1. Wave Maker and Numerical Dissipation 


Many of the issues with boundary conditions came as a result of trying to have a 
boundary that would produce a wave field, and also have a constant current flow. This 
study facilitated the creation of an internal wave maker. The wave maker is rather simple 
in principle in that it requires an amplitude, period, and depth input parameter. It then 
creates a wave field from within the tank boundaries with propagation in both the positive 
and negative x direction. By having a wave source contained within the boundary, the 
issue of dual dynamic conditions at the east and west boundaries was solved. However, 
the simplistic wave maker did not produce waves at the exact amplitude as the input, 


therefore adjustments to initial conditions were needed. 


Initial tests of the model with waves only in the tank were conducted with a 0.1 m 
resolution in the x and y directions. The wave maker in the initial test introduced a wave 
field of 0.03 m amplitude and 1.2 s period. At this resolution, the internal linear wave 
maker is unable to produce a constant amplitude wave. Issues of numerical dissipation, 
seen in model data, caused amplitude damping that does not support wave theory (Figure 
2). Amplitude damping occurred within 10 m of the wave maker in both a square tank, 


and the test tank with the false wall simulating an inlet at this resolution. 
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Figure 2. Wave height distribution through the tank at low resolution (0.1m) a.) 
Square tank. b.) Tank with false inlet wall. Vertical lines indicate 
first narrowing wall of the tank domain. 


By increasing the resolution to 0.025 in the x direction and 0.05 in the y-direction 
the wave maker was able to produce steady waves that maintained amplitude with 
minimal damping through a square tank as seen in the wave height plotted through the 
model tank (Figure 3, blue line). However, amplitude damping was still present in the 
tank with the inlet wall (red line) as a result of a boundary condition issue to be described 
hereafter. Each case shows a consistent wave through the wide area of the tank; however, 


the red profile begins to dampen considerably upon reaching the inlet wall at 39.4 m. 
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Figure 3. | Wave height at high resolution for waves in a square tank 
with no current (blue), and waves in inlet tank, no current (red). 


Numerical dissipation is caused by a numerical approximation of the partial 
differential equations assuming the theoretical equations are exact. For example, using a 
first order finite differencing scheme will truncate equations at the second order (dx’), 
The truncated terms may be diffusive or dispersive. In this model, the terms are diffusive. 
NHWAVE uses a numerical approximation to the partial differential equations governing 
the internal wave maker in the tank. This process causes a dissipation of the wave down 
the tank that is a result of the numerics of the model rather than an actual physical 
process. In the model results, the wave maker showed amplitude damping throughout the 


tank that improved with increasing the horizontal resolution. 
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yas Boundary Conditions 


In order to model the CKO2 tank, great care was taken in developing the 
NHWAVE boundary conditions as mentioned previously. Having an internal wave maker 
was beneficial in simplifying the boundary inputs, however, other issues arose. Because 
the waves were created inside the domain, waves propagate in both directions and need to 
be dissipated at the boundaries to negate reflection from the east and west walls. 
Absorption of mass flux and energy at the boundary is also a concern for the current 
imposed in the domain. If the current is allowed to reflect off the boundaries, oscillations 


in the tank will form and create issues in the data. 


Sponge layers on either side of the tank were created to dissipate wave energy and 
mass flux from the current. Throughout the modeling process modifications to the source 
code, to include dimension and decay coefficients, were needed in order to negate as 
much reflection from the tank “walls” as possible. By creating a low-pass filter imbedded 
in the sponge layers, oscillations and reflections were decreased to almost zero. The filter 
allows the current to pass through at a constant speed while filtering out any wave action. 
Wave energy is absorbed into the sponge layers at the east and west boundaries (Figure 
4). Development of this aspect of the model is crucial to future iterations and uses of the 


model in experimental tank applications. 
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Figure 4. _ Low pass filtered sponge layers showing wave dissipation in the 
sponge layer with a .03 m 1.2 s wave, no reflection. 
Sponge layers are? 7 m wide. 


No-slip boundaries were used on the north and south walls of the tank. This was 
done to reduce any change in current velocities or the wave field in the cross-tank 
dimension. Also, the bottom was defined as frictionless to maintain the constant vertical 
distribution of current assumption. The vertical and cross-tank distribution of current for 
the test cases shows little variability (Figure 5). This allows one-dimensional theories to 


be applied to the CK02 and model results. 
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Figure 5. Vertical current distribution at a) 31.5 mand b) 33 m.c) 
Horizontal distribution across the tank at location 25 m. 


3. Computational Efficiency 


NHWAVE is set up to run in either parallel mode, using the message-passing 
interface (MPI), or in serial mode where a single node of computing is used. For this 
study, the high-powered computers in the Mills cluster at the University of Delaware 
were accessed to run the model. The model was programmed in Linux and Fortran. Due 
to nodal structure, computational efficiency was greatly increased by running in parallel 
mode. This was especially evident as grid resolution in the model was increased. Due to 
the nodal structure of the computers, computational hours were shared with the Coastal 
Engineering group at UDel. At times, finding available nodes to run model simulations 


was difficult and caused delays in data acquisition. 
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To compound the issues of nodal effects in high-powered computing, the Mills 
cluster suffered a catastrophic failure. Of the six hard disks comprising the cluster 
system, data can be preserved if two or fewer disks suffer failure. In the process of 
running and compiling this model, three hard disks were destroyed, this caused a 
catastrophic failure of the supercomputing capability at UDel and data along with model 


updates were lost. Recreating the model data was both difficult and time consuming. 
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IV. EXPERIMENTAL APPROACH 


A. INTRODUCTION 


Measurements from the accompanying data report to CK02 are used to test 
NHWAVE model output for cases of refractive shoaling through a tank of varying width. 
Each simulation was done using monochromatic waves on a constant opposing tidal flux. 
A two-dimensional wave tank domain modeled after the CK02 tank experiment was 
created to simulate a river inlet. The goal was to see how well NHWAVE predicts wave- 
current dynamics in regions where wave-blocking occurs and compare results to those of 


CK02 in the hope of verifying the application of NHWAVE to wave blocking cases. 


B. TEST TANK CONFIGURATION 


The CK02 data set was acquired in a 30.0 m long recirculating wave flume at the 
Center for Applied Coastal Research at the University of Delaware. The width of the 
flume is 0.60 m wide with a constant depth of 0.50 m (Figure 6). The coordinate system 
is oriented so that the positive horizontal axis, x, is the direction of wave propagation, 
x=0 at the location of wave generation. The y-coordinate axis has an origin on the 
continuous sidewall (right side) and points toward the narrowed wall that was created to 


represent the inlet. The z-axis is positive upwards with z=0 at the still water level. 


A false wall was implemented starting at x=12.4 m to gradually narrow the tank 
to 0.36 m by 15.2 m in the tank. The narrowing of the channel was designed to represent 
a tidal inlet where energy density and current velocity would increase in the narrow 
section of the tank. This design, coupled with a constant opposing current, was intended 
to block waves with periods up to 1.3 s in a 0.5 m depth. The depth-averaged current 
throughout the flume was increased through the narrow shoal region, but maintained 
relative depth uniformity. After 4.9 m of the narrow section, the tank width gradually 
widens at five degrees to the original dimensions negating flow separation at the 


transition. 
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Figure 6. Schematic for Chawla and Kirby tank experiment 
(After Chawla and Kirby 2002). Not to Scale. 


Waves were generated on the west side, (x=0). The current was created by 
drawing water from behind the perforated wave maker and pumping it to the opposite 
(Eastern) end of the tank. Care was taken to ensure there were no large-scale eddies. A 
complete description of the mechanics of the wave tank and the experimental method in 


the lab are found in Chawla and Kirby (1999). 


For the CK02 data, capacitance wave gauges were used to measure the free 
surface elevation and Acoustic Doppler Velocimeters (ADVs) measured fluid velocities. 
In the model applications, free surface elevation and fluid velocities were calculated at 


approximately the same locations as the CK02 instruments (Figure 6). 


C. MODEL SET-UP 
1. Domain 


NHWAVE was initialized with a domain that was altered from the original tank 
experiments to account for model boundary condition constraints. The bottom grid was 
produced in a Matlab environment to simulate a tank that measured 70 m X 0.6 m (Figure 
7). By increasing the length of the tank, effects from the wave-maker were reduced. The 
coordinate axis are the same as in the CK02 experiment. The computational grid was 


configured with 0.025 m resolution in the x-direction and 0.05 m resolution in the y- 
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direction. The depth was set at 0.5 m and the model was run with five vertical layers. A 
false wall was initiated at 39.4 m increasing in width from the tank sidewall boundary to 
narrow the channel to 0.4 m at 42.2 m in the domain. The simulated inlet then extends for 
10 m before expanding at the same ratio to the original with of the tank, see Figure 3. 
A time step of 0.1 seconds was used for computational efficiency while maintaining high 
enough temporal resolution. Output was computed at gauge locations closely resembling 


the locations in CK02. 


Approximate Gauge Locations —_ 


COKODOOD OOD OOD OOD OOD0OROO00000 


Tidal flux 


Wi 





Figure 7. | Schematic for NHWAVE model setup. Not to scale. 


The artificial inlet wall was moved to many different locations in the tank and 
changed in dimension to attempt to replicate the currents and wave amplitudes present in 
CK02. The wall was defined with a no-slip boundary, which should have produced wave 
reflection and amplification in the tank when a current was not included. However, as the 
model results indicate, there was significant amplitude damping through the narrowed 
inlet (Figure 3). After increasing resolution to account for numerical dissipation, it was 
evident that the boundary conditions were not resolved for the model domain when a 


false wall was implemented. 


Previous versions of NHWAVE were run in a domain replicating bathymetry 
present in the field. The boundaries in these cases are dissipative in that landmasses are 
absorb wave energy. This land masking allows the model to run with dissipation effects 


that would be expected in the environment. In the tank domain the depth of the tank was 
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0.5m. The false inlet wall was defined by a trapezoid that stood 1.0 m above the still 
water depth. This involves the model needing to account for wetting and drying of the 
boundary. The non-hydrostatic approximation allowed for pressure perturbations and 
vertical velocities to be introduced in the domain. In the tank model runs, the pressure 
field showed errors in the presence of dry points in the domain. This boundary issue will 


need to be resolved in the future development of the model. 


2s Wave Maker 


Waves were created by a linear wave maker located 15 m in the tank. Due to the 
simplicity of the wave-maker, model input was adjusted to replicate wave amplitude at 
the first gauge in the experimental data. In CK02 the initial conditions were defined at the 
first gauge location that varied between 4.6 m to 5.8 m before the narrowest part of the 
channel; the initial conditions in the model were defined at similar distances from the 


narrowed inlet but implemented to replicate CK02 experimental data. 


3. Sponge Layers 


Sponge layers 7 m thick were placed at each end. The sponge layers were made 
thick and with a low decay rate to allow waves and currents to reach the boundary with 
little reflection back into the tank from the east/west boundary. This tank set-up with 
appropriate sponge layers allows simulating a refraction only inlet where wave blocking 


occurs close to the narrow part of the “inlet.” 


4, Current Flow 


The mean current was assumed uniform in the vertical and cross-shore, due to the 


no-slip boundary verified previously. The current was given by 
=— g (16) 


where b is the width of the channel and Q is the volume flux prescribed as a tidal flux in 
the input boundary conditions. In this equation, boundary layer effects are ignored so that 
the value of Q was adjusted to 0.05 m°/s to match the data. The measured and calculated 


horizontal current profile through the narrowing of the tank is presented in (Figure 8). 
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The profile agrees well with theory in the transition to the inlet and the faster laminar 
flow located there. There is an increase in the current speed just prior to the inlet opening 
back to the width of the channel. This jump could lead to inconsistencies in data at the 
inlet transition area. This deviation from the tank data is due to a pressure field imbalance 


at the inlet boundary to be explained in Chapter V. 
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Figure 8. Horizontal current distribution through the narrowing of the inlet. 
(Model output circles (0), and theory solid line.) 


As determined in the development of the model, vertical and cross-tank horizontal 
velocity mean profiles through 400 s in five different locations within the tank showed 
little variation, and thus the uniform assumption was valid. This allows for application of 
equations in a one-dimensional sense. Velocity profiles throughout the tank show relative 


depth uniformity with an average current of 0.2 m/s outside the inlet and 0.32 m/s 


va 


throughout the narrowed section (Figure 9). Decreasing the current from CK02 


observations allowed a slower transition through the inlet and less instability in the model 
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Figure 9. Vertical current distribution in five locations along the tank. 


The current in the model tank showed inconsistencies through the narrowed inlet. 
Laminar flow was present in the wide area of the tank at both the east and west side just 
after the sponge layers. As the current flows over the narrowing section of the wall from 
the east, there is marked instability. Laminar flow is then re-established through the 
narrow section, this is where model conditions were measured from for data comparison 
(Figure 10). The surface current becomes faster than the current at the bottom of the tank. 
This creates a bore-like situation where one shouldn’t be present and suggest that there is 


a boundary problem. The current then “jumps’ before exiting the inlet. This could be the 
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cause for what appears to be wave blocking in the model data, however it is a function of 


problems with the pressure field around the false wall. 


In examining the current only case the surface also shows shape change over the 
opening of the inlet (Figure 10). The inconsistencies over the eastern narrowing are 
present, as well as a dip in elevation at the opening of the inlet. Froude numbers in the 
narrowed inlet were 0.14 and decreased to 0.07 in the less constrained flow. These 
numbers do not indicate an approach to critical flow and therefore, super-critical flow 
should be excluded from consideration for the formation of a bore-like element. These 
inconsistencies in the data perpetuated more development of the numerical scheme for 


boundary conditions. 
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Figure 10. Surface elevation in the current only case (green). 
Current profiles at surface (red), and at tank bottom (blue). 


The wet/dry boundary issue also causes vertical velocities to appear with changes 


in the pressure field. The laminar flow should show no vertical accelerations whether in 
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the inlet or in the wide area of the tank. As the current flows through the inlet, the orbital 
velocities show vertical flow, an unexpected result in the model (Figure 11). This occurs 
both at the eastern side of the inlet, where flow is restricted, and also at the opening of the 
inlet to dissipative laminar flow. This result is shown through all layers of the current 


flow. 
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Figure 11. Vertical velocity instabilities in the current through the inlet. Multiple colors 
indicate multiple layers from magenta at the surface, to blue, red, and then green 
at the bottom. 
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V. WAVE BLOCKING TESTS 


A. INTRODUCTION 


The extensive data that were collected for CK02 is compared to NHWAVE model 
data. In the CK02 test cases, wave periods ranged from 1.2 to 1.6 s. The waves created in 
the experiment were intermediate-depth water waves with kh values outside the narrow 
channel ranging from 1.35 to 2.4, as determined using a Doppler-shifted linear dispersion 
relation (Chawla and Kirby 2002). The kh values increase inside the narrow channel 
where currents are stronger. By narrowing the channel in the experiment, waves were 
shown to shoal with or without an opposing current. Conditions meant to replicate this 
data were explored in the model. The results of the NHWAVE model runs were not 
accurate in predicting the wave blocking and breaking conditions as seen in CK02. 
Constraints in modeling in a tank domain have been found to cause instabilities and false 


results in the model. 


Each case for monochromatic waves was divided into two parts, small amplitude 
where waves should reflect from the blocking point without breaking, and larger 
amplitude waves where breaking occurs. The study of the amplitude envelope through 
the blocking region as well as the energy transformation and dissipation through the 
region was accomplished by implementing a range of wave heights with different wave 
periods in the experimental and model initial conditions. Each model run was conducted 
to attempt to mimic the wave blocking results in CK02, however, the model was 


unsuccessful in resolving the wave components with the current numerical scheme. 


B. THEORETICAL EXPECTATIONS AND CK02 RESULTS 
1. Monochromatic Small Amplitude Reflected and Blocked Wave Tests 


In monochromatic, small amplitude tests, wave energy cannot propagate beyond 
the blocking point and no energy is lost due to breaking at or before the blocking point; 
therefore, there is wave reflection. From previous studies, these reflected waves are 
unique in that the phase speed of the waves continues to move against the current, 


however, the group velocity moves with the current with the wave energy (Chawla 1999). 
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For the monochromatic wave blocking tests, three different periods were used with 
increasing wave amplitude to study the wave amplitude dynamics through the blocking 


region and were compared to CK02 data. 


In linear wave theory, if there is no wave breaking in the domain and the wave 
amplitude and steepness is small enough that waves will be reflected from the blocking 
point with a different wave number. The amplitude envelope of these waves through the 
blocking point resembles an Airy function (Peregrine 1976). CKO02 data verified this 


results for small amplitude monochromatic wave tests. 


In the experimental data, the change in the amplitude envelope is clearly visible in 
the smallest wave cases. The test cases presented here are for the 1.2 s period. Wave 
height measurements are made in relation to the blocking point with x=0 indicating the 
narrowest part of the test tank. As the amplitude increased and approached the breaking 
steepness criteria as defined by Miche (1951), the node and antinodal resemblance to an 
Airy function is diminished (Figure 12). Non-linear effects change the steepness of the 
wave reflection response, which can be seen by the increasing deviation in the primary 


wave’s peak and slope as the amplitude increases. 
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Figure 12. Comparison between measured CKO2 data and Airy function 
for small amplitude monochromatic waves. Amplitude increases 
from 1.3 cm in test 2 to 1.6 cm in test 5 (From Chawla 1999). 


2s Monochromatic Breaking Wave Tests 


In CK02, each test was conducted similar to the small amplitude conditions in that 
multiple amplitudes were examined for different wave periods. A simple numerical 
model was created to study the wave breaking in this tank configuration. The bore model 
using third order Stokes relations were able to handle the location of the breaking and 
blocking of the wave field well. However, it was insufficient to predict the side-band 
instabilities and energy shifts seen in the data. The hope was that NHWAVE would be 
able to predict these. 
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a. Energy Spectra 


One of the main findings in CK02 was the frequency downshift in the 
energy spectrum of monochromatic breaking waves. This shift is a result of the nonlinear 
side band instabilities. These are dependent on the frequency and amplitude of the wave 
field and effect the wave propagation greatly (Benjamin and Feir 1967). As the waves 
become steeper on the opposing currents and approach the blocking point, the group 
velocity approaches zero and wave energy travels very slowly (Lai et al. 1989). Therefore 
energy can be shifted to a lower sideband, which requires a higher blocking current. This 
will allow the wave components in this lower side band to travel into the narrow inlet 
while the higher frequency waves are blocked, and further increase the energy in the 
lower side-band as a function of distance through the inlet. The result is that the primary 
wave component has reached its blocking condition, while the lower sideband propagates 


through the channel. 


The spectra from CK02 wave breaking conditions show the frequency 
shift to lower sidebands (Figure 13). For probe locations around the blocking point, 
indicated in CK02 as x=0, the primary (P), upper (U), and lower (L) sidebands show 
modification of the frequency spectra. Energy in the upper and primary sideband 
decreases with progression into the inlet while there is marked lower sideband growth. 
This supports the theory of higher frequency energy being blocked while lower sideband 
waves are allowed to propagate through the inlet. The non-linear aspect of the lower 


sideband growth was not predicted by the simplified model used in CK02. 
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Figure 13. Frequency spectra for a wave breaking test from CK02 with period=1.26 s 
and amplitude=0.126m (L=lower side band, P=primary side band, 
U=upper sideband) (From Chawla 1999). 


C. NHWAVE PRELIMINARY RESULTS 


Model output of surface height, n, and velocities in the u, v, and w directions were 
output in the model over the domain and through the five vertical layers in the simulated 
tank. Comparison tests for CK02 data were unsuccessful. However, recent developments 


have made some wave blocking cases in a tank with a narrowed inlet possible. 


1. NHWAVE Performance 


Each case of monochromatic wave blocking and wave breaking was initially 
attempted in NHWAVE for comparison to data from CK02. Wave height distributions 
through the blocking point show vast differences from theory and experimental data due 
to the unresolved boundary conditions. In some cases NHWAVE showed a wave height 
distribution that appeared to have wave blocking for the range of small amplitude waves 
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tested. In these cases, the reflected wave was not apparent in any of the measurements. 
Coupled with the wave dissipation caused by the pressure field inaccuracy through the 


inlet, these results do not verify the model. 


Larger amplitude waves cases were also examined using the wave height 
distribution to determine if the waves broke at, or just prior to, the blocking point of the 
inlet. In those tests where a steep drop in amplitude around the blocking point indicated 
that waves were breaking, wave period and spectra were examined and compared to 
CKO02 data. The model shows significant deviation from the expected wave height 


distribution and energy spectra through the blocking point and into the inlet. 


2s Recent Model Developments 


Preliminary results from a new iteration of NHWAVE are promising for wave 
blocking cases. For these cases a new tank domain that measured 35 m long by 0.6 m 
wide was created (Figure 14). The tank has a simulated inlet 7.5 m long centered at 20 m 
in the tank. The inlet narrowed the overall width of the tank to 0.3 m. For the new test 
cases, resolution was 0.025 in x and 0.05 in y. This was shown previously to negate 
numerical dissipation issues from low resolution test runs. Sponge layers were 5 m thick 
on either end of the tank and had the same decay coefficients and low pass filter as used 


previously. 
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Figure 14. Simplified test tank for preliminary model results. (Not to Scale) 
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By implementing a new Poisson solver for the issue of the land masking and the 
wet/dry boundary, the pressure field inconsistencies have been resolved. In a simple tank 
with an inlet wall halfway down the “tank” wave amplitude increases through the inlet 
as would be expected before decreasing due to diffraction as the tank widens again 


(Figure 15). 
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Figure 15. Wave height in simplified tank. Inlet wall location is indicated by the 
trapezoid, the wave maker is located at the vertical line. 


a. Preliminary Wave Blocking Results 


For small amplitude waves, blocking conditions were met in the 
NHWAVE data. Waves shoaled slightly as they approached the wall and decrease in 
amplitude due to the opposing current at the end of the inlet. As the current is increased 
from 0.28 m/s (Figure 16a) in the tank to 0.35 m/s (Figure 16b), the blocking criteria is 


reached for the test wave of period 1.3 s and initial amplitude of 0.09 m. With increased 


a7 


current, the shoaling is more pronounced and the wave steepness is greater. This causes 
the wave to break as seen in the rapid decrease in wave amplitude. Therefore, for small 
amplitude waves, the effect of increased current through the inlet shows blocking criteria 


being met for the given wave parameters. 
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Figure 16. Small amplitude wave blocking and reflection. Solid lines indicate the 
boundaries of the inlet wall. a) Opposing current=0.28 m/s b) Opposing 
current=0.35 m/s. Vertical lines indicate the boundaries of the inlet wall. 


Similarly, in a large amplitude cases, with a period of 1.3 s and amplitude 
of 30 cm, waves show an increase in amplitude just prior to the narrowing of the inlet and 
breaking at the narrowest section for the larger current. Again implementing an opposing 
current of 0.28m/s caused the waves to shoal against an opposing current, but not break 
(Figure 17a). With increased speed, 0.35 m/s, the waves broke upon entering the 
narrowed inlet (Figure 17b). This is consistent with the theory presented here as well as 


the results from CK02. 
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Figure 17. Wave height distribution for a large amplitude wave breaking case. 
Vertical lines indicate the boundaries of the inlet wall. 


Time series analysis for energy spectra were not available due to the 
shortened data record for the simplified cases. Further statistical analysis is needed to 
determine if NHWAVE can be validated in comparison to CK02 for wave blocking 


Cases. 


3. Extension to Random Wave Spectra 


The validation attempts for NHWAVE were conducted with monochromatic 
waves only. One of the benefits of the CK02 data was examining the changes in the 
frequency spectra through the breaking and blocking region with wave groups. Extension 
to more realistic spectra for a better understanding of the environment is needed in 
NHWAVE. Currently, efforts to include implementing Jonswap spectra are being 


attempted. However, results are not reliable at this time. 
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VI. CONCLUSIONS 


The complex dynamics of wave blocking are difficult to model. Current wave 
averaged models under predict the non-linear processes that occur within short distances 
of a wave blocking location. Two models were recently developed to model fully 
dispersive wave processes and were explored for their use in river inlets where blocking 
conditions exist. The testing of both SWASH and NHWAVE proved boundary conditions 
involved in a tank setting necessitated modifications to the model code and numerical 
processes therein. Multiple iterations of NHWAVE were developed for this application 


and resolution of boundary condition issues proved difficult, but in the end, successful. 


Model testing using a series of monochromatic wave cases on opposing blocking 
currents was conducted with NHWAVE. Two classifications of tests were presented, 
those with small amplitudes where waves were reflected, and those with higher 
amplitudes that demonstrated wave breaking prior to the blocking point. These tests were 


compared to data from CK02. 


Data from CK02 were able to show monochromatic wave blocking and wave 
breaking on opposing currents. Waves that were reflected at the blocking point showed a 
distribution close to an Airy function for small amplitude waves. This is a close 
approximation to linear theory that breaks down as wave amplitude increases. In larger 
amplitude wave breaking cases non-linear effects are increased and seen in the energy 
spectra. The changes in the energy spectra indicate unstable sideband growth that is 
difficult to model. NHWAVE in its current state is unable to replicate the CK02 tank 
experiment due to boundary conditions that cause amplitude dissipation through the 


narrowing of the inlet. 


Future work on random wave field input to NHWAVE as a Jonswap spectrum is 
needed. At the moment, the model is unable to handle a random wave field or wave 


group that would be more pertinent to field data. 


NHWAVE, while still under development, allows us to explore many conditions 


that are not always possible with field or lab measurements and shows promise as an 
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accurate predictor of wave-current interactions and wave blocking. In comparison to 
CK02 lab data, there are shortfalls in the numerics of the model that are greatly affected 
by boundary conditions. Preliminary results with improved equations for boundary 
conditions show promise for future implementation of the model. More development will 
lead to improved understanding of the wave blocking conditions in river inlets. The 
future implementation of a non-hydrostatic model for the intense kinematic conditions in 
wave-cutrent interactions and wave blocking scenarios will be a valuable addition to our 


knowledge base of these turbulent areas. 
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